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Abstract 

This  thesis  addresses  the  prob].em  of  estimating  the 
target  angle  with  respect  to  the  boresite  of  an  Amplitude- 
Comparison  Monopulse  Radar.  The  maximum  likelihood  esti¬ 
mate  is  used  to  track  the  target  and  produce  a  mean-square 
error.  This  error  is  approximated  from  a  computer  simu¬ 
lated  tracking  loop  and  then  compared  to  the  Cramer-Rao 
and  Ziv-Zakai  bounds. 

At  high  signal-to~noise  ratios  (SNR) ,  the  Cramer- 
Rao  bound  is  useful  in  lower  bounding  the  mean-square  error. 
At  low  SNR  levels  a  phenomena  known  as  the  threshold  effect 
occurs  and  the  Cramer-Rao  bound  becomes  unreliable  as  a 
lower  bound.  The  Ziv-Zakai  bound  is  a  tighter  bound  and 
can  be  used  to  lower  bound  the  mean-square  errcr  at  very 
low  SNR  levels.  It  also  proves  useful  in  determining  the 
SNR  level  at  which  the  threshold  effect  occurs. 
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LOWER  BOUNDS  APPLIED  TO  THE  MEAN-SQUARE  TRACKING 


ERROR  OF  AN  AMPLITUDE-COMPARISON 
MONOPULSE  RADAR 


I.  Introduction 


Background 

A  monopulse  radar  tracks  by  maintaining  the  radar 
boresite  on  the  target.  As  the  target  moves  off  of  bore- 
site  an  error  voltage  is  developed  proportional  to  the 
angle  off  boresite.  This  information  is  fed  to  the  rad  ir 
angle  servo  which  in  turn  steers  the  boresite  toward  the 
target.  Theoretically,  this  process  is  done  using  infor¬ 
mation  gained  from  a  single  pulse.  Breaklock  occurs  at  a 
certain  signal-to-noise  ratio  (SNR)  per  pulse  when  the 
information  gained  from  the  pulse  is  not  sufficient  to 
steer  the  boresite  to  the  target  accurately.  Consequently, 
the  target  moves  out  of  the  antenna  beam  width. 

At  a  certain  minimum  SNR  level,  a  phenomena  known 
as  the  threshold  effect  occurs.  Consider  the  example  of  a 
probability  density  function  (pdf)  of  an  observation  z. 
Assume  that  the  density  function  has  several  maximums,  the 
largest  of  which  is  at  the  center.  The  maximum  likelihood 
(ML)  estimate  of  z  is  that  value  of  z  which  maximizes  the 
pdf  (Ref  6] . 
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For  the  case  of  no  noise  in  the  observation  z, 


the  maximum  is  at  the  origin  (Figure  1) .  For  the  low  noise 
situation,  the  maximum  is  shiftec'  a  small  distance  from  the 
origin  which  is  proportional  to  the  standard  deviation 
of  the  additive  noise.  When  the  ncise  becomes  comparable 
to  the  mean  of  the  pdf,  not  only  is  the  local  maximum  near 
the  origin  shifted  but  it  is  highly  probable  that  some 
other  local  maximum  of  the  pdf  becomes  the  largest.  Now 
the  error  is  much  larger  than  the  error  caused  by  the  shift 
ing  of  the  central  local  maximum.  These  large  errors  are 
referred  to  as  anomalous  errors.  When  anomalous  errors 
occur,  the  system  behaves  nonlinearly  and  its  performance 
degrades  rapidly.  This  is  what  is  known  as  the  threshold 
effect  [Ref  11] . 


Fig.  1.  Probability  Density  Function  of  Observation  z 


The  incoming  target  return  is  processed  in  the 
monopulse  radar  by  weighting  the  amplitude  of  the  return  by 
the  antenna  gain  factor  of  its  sum  and  difference  beams. 
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This  amplitude  weighting  is  proportional  to  the  target 
angle  off  of  boresite  and  is  used  to  steer  the  radar  bore- 
site  to  the  target.  Because  of  random  channel  noise  how¬ 
ever,  the  boresite  will  be  steered  to  the  vicinity  of  the 
actual  target  position.  The  vicinity  is  described  in  terms 
of  the  mean- square  error  (MSE) .  It  will  be  shown  hereafter 
that  as  the  SNR  level  decreases,  the  MSE  increases. 

The  threshold  can  be  quantified  by  describing  a 
mean-square  tracking  error,  driving  its  lower  bound  and 
then  comparing  the  two.  It  can  be  shown  [Ref  2]  that  the 
Cramer-Rao  (CR)  bound  for  the  variance  of  the  probability 
density  of  the  angle  off  boresite  estimate  is  proportional 
to  the  beam-width  and  SNR. 


Var ( 0 ) 


BW 
“  SNR 


(1) 


This  bound  tells  us  that  no  matter  what  SNR,  the 
variance  is  kept  small  by  maintaining  a  narrow  beam-width. 
However,  the  CR  bound  is  only  useful  when  the  SNR  level 
is  well  above  the  system  threshold.  Below  threshold  the 
CR  bound  becomes  inaccurate  and  proves  unable  to  predict 
breaklock.  A  much  tighter  bound  (one  good  at  low  SNR)  is 
needed  to  quantify  and  predict  the  system  threshold  or 
breaklock  point. 


Problem  and  Scope 


The  purpose  of  this  thesis  is  to  develop  a  maximum 
likelihood  estimate  of  the  target  angle  off  boresite  and 
use  it  to  track  the  target.  From  this  procedure  the  MSE 
can  be  derived  through  simulation  and  used  to  determine 
at  what  SNR  the  MSE  becomes  large  enough  to  be  considered 
a  breaklock  point.  The  simulation  results  will  then  be 
compared  to  the  CR  and  Ziv-Zakai  bounds  to  see  which  bound 
will  provide  an  accurate  prediction  of  the  breaklock  point. 

The  Ziv-Zakai  (ZZ)  bound  is  designed  to  be  tighter 
at  low  SNR  levels  [Ref  11] .  The  bound  was  derived  by  con¬ 
sidering  the  suboptimal  detection  procedure  which  involves 
an  estimate  of  one  of  two  possible  values  of  a  parameter. 

A  criteria  was  developed  for  determining  which  value  of 
the  parameter  was  estimated.  This  approach  is  called  the 
Estimation  Theorist's  Approach.  The  probability  of  error 
associated  with  this  approach  is  lower  bounded  by  the  proba¬ 
bility  of  error  achieved  by  the  optimal  detection  scheme 
(which  is  referred  to  as  the  Hypothesis  Tester's  Approach) . 
This  entire  process  resulted  in  new  bounds  based  on  known 
results  in  detection  theory.  The  advantage  of  the  ZZ  bound 
is  that  it  is  independent  of  biases  associated  with  esti¬ 
mators  and  is  therefore  tighter  than  other  bounds  [Ref  11] . 

The  ZZ  bound  was  modified  to  fit  a  problem  involv¬ 
ing  the  bearing  estimate  of  a  linear  array  [Ref  8],  This 
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modified  version  is  adapted  and  used  to  derive  a  lower 
bound  for  the  MSE  for  this  problem. 

Assumptions 

This  problem  is  approached  assuming  a  one¬ 
dimensional  Amplitude-Comparison  Monopulse  Radar  (i.e., 
azimuth  case  only) .  The  channel  noise  is  assumed  to  be 
white  gaussian  noise  (WGN)  independent  from  channel  to 
channel  and  pulse  to  pulse  with  zero  mean  and  variance 

N  /2.  The  antenna  functions  used  to  describe  the  sum  and 
o 

difference  beams  are  derived  assuming  <.he  two  radar 
antennas  are  point  sources  in  a  linear  array.  The  ampli¬ 
tude  of  the  return  signal,  A,  is  assumed  to  be  random  and 
unknown  since  the  target  radar  cross-section  is  unknown. 
However,  A  can  be  estimated  as  will  be  shown  and  used  in 
the  estimation  of  the  target  angle  relative  to  boresite. 

The  radar's  3  "3  beam-width  is  3°  and  the  PRF  is 

200  pps. 

A  stationary  target  will  be  tracked  first  after 
which  a  moving  target  will  be  considered. 

General  Approach 

The  first  step  is  to  develop  adequate  antenna  func¬ 
tions  which  can  weight  the  return  signal  amplitude  propor¬ 
tional  to  the  target  angle  relative  to  boresite.  The 
antenna  functions  are  used  to  form  the  observation  r(t). 

The  product  of  the  return  signal  s(t),  its  amplitude  A, 
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and  the  antenna  functic  i,  added  to  the  channel  noise  makes 
up  an  observation. 

r(t)  =  As(t)g(0)  +  n(t) 

It  is  shown  that  the  one-dimensional  monopulse 
radar  has  two  channels  and  will  therefore  have  two  observa¬ 
tions  r^(t)  and  r2(t).  With  these  observations  and  their 
joint  pdf,  the  maximum  likelihood  estimate  can  be  derived 
by  finding 


£gjp  p(r1r2/A,6)  =  0  (2) 

where  p(r^r2/A,0)  is  the  joint  pdf  of  r^  and  r2»  A  is  the 
amplitude  of  the  return  signal  and  0  is  the  target  angle 
relative  to  boresite. 

This  results  in  an  equation  that  is  a  function  of 

A 

0,  the  actual  target  position  relative  to  boresite,  and  0, 
the  estimated  target  position.  It  is  then  shown  that  this 

A 

equation  is  zero  when  0  =0,  or  when  the  radar  boresite 
is  pointing  at  the  actual  target  position.  This  equation 
is  called  the  error  equation  and  is  used  to  steer  the 
antenna  toward  the  target. 

Once  the  error  equation  is  obtained,  the  ZZ  bound 
is  derived  following  the  steps  outlined  in  reference  [8] . 

At  this  point,  a  conputer  simulated  tracking  loop 
is  designed  and  implemented.  The  tracking  loop  will  yield 
an  error  between  the  actual  target  position  and  the 
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boresite  position.  This  error  is  developed  into  an  MSE 
and  then  plotted  for  different  SNR  levels.  The  MSE  will 
be  a  function  of  time  since  an  error  will  develop  for  each 
pulse  or  observation  received.  The  breaklock  point  will 
occur  when  the  MSE  grows  unbounded  as  the  number  of  pulses 
(or  observations)  increases.  The  level  of  SNR  where  this 
occurs  is  compared  to  the  Cramer-Rao  and  Ziv-Zakai  bounds 
to  quantify  the  system  threshold  or  breaklock  point. 

Both  the  stationary  and  nonstationary  targets  are 
used  to  determine  breaklock. 
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II.  The  Radar  Model 


The  Antenna  Functions 

A  simple  block  diagram  of  a  one-dimensional  mono¬ 
pulse  radar  is  shown  in  Figure  2. 


Fig.  2.  Simple  Block  Diagram  of  Radar 

The  patterns  of  antennas  A  and  B  are  combined  so 
as  to  produce  two  beams.  The  sum  beam,  A+B,  is  represented 
as  an  antenna  function  g^,(9).  Likewise,  the  difference 
beam,  A-B,  is  represented  as  g^(0).  The  antenna  functions 
are  used  to  describe  the  amplitude  weighting  of  the  target 
return.  The  maximum  gain  occurs  at  the  antenna  boresite 
(0=0)  and  decreases  as  0  increases. 

To  derive  the  antenna  functions,  the  two  antennas 
are  considered  as  point  sources  in  a  two-element  linear 
array  (see  Figure  3). 


Fig.  3.  Two  Point  Sources  in  the  x-z  Plane 


A  plane  wave  approaching  from  an  arbitrary  direc¬ 
tion  in  the  x-2  plane  can  be  described  as 

(2) 
(3) 

where  Eq  is  the  amplitude  of  the  plane  wave  and 


E.  =  E~e 
A  o 


eb  *  V 


jkB-r 


K.  ss  k  X  X 

A  o  o 


KB  “  'Vox 


and 


r  =  x  sin  0  cos  <p  +  y  sin  0  sin  <p  +  z  cos  0 


Performing  the  dot  product  and  substituting  back  into 
equation^  (2)  and  (3)  gives 


E.  =  E  e 
A  o 


jk,x  sin  0 
o  o 


(4) 


(5) 


EB  =  EoG 


-jk  x  sin  0 
o  o 


Equations  (4)  and  (5)  can  now  be  used  to  form  the  sum 
and  difference  beams  of  the  monopulse  radar. 


jk  x  sin  0  -jk  x  sin  0 
ea  +  Efi  =  EQ(e  °  °  +  e  °  °  )  (6) 

jk  x  sin  0  -jk  x  sin  0 
EA  -  EB  -  Eo(e  °  °  -  e  0  °  )  (7) 


Using  Euler's  equation  and  taking  the  magnitude  yields 
[Ref  8] 


!ea 

+  ebI  " 

2Eo 

cos  u 

(8) 

Iea 

-  ebI  ■ 

2Eo 

sin  u 

(9) 

where  u  =  k.  x 

n  A 

sin  0 . 

o  o 


The  antenna  functions  for  the  radar  model  can  be 
obtained  by  normalizing  equations  (8)  and  (9)  to  produce 


g  (u)  =  cos  u 

Lt 

(10) 

g.  (u)  =  sin  u 

A 

(ID 

These  functions  represent  the  monopulse  sum  and 
difference  beams.  The  patterns  from  the  two  antennas  are 
approximated  in  Figure  4.  To  obtain  the  sum  beam  these 
two  patterns  are  summed  and  yield  the  pattern  shown  in 
Figure  5a.  The  difference  beam  is  obtained  by  subtracting 
the  two  patterns  which  yields  the  pattern  in  Figure  5b. 
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j.f  the  boresite  is  pointing  directly  at  the  target 
g£ (u)  =  cos  (sin(O) )  =  I 

and 

g^(u)  =  sin  (sin(O))  =  0 

Now  that  the  antenna  functions  have  been  estab¬ 
lished,  a  model  for  estimating  the  angle  off  boresite  can 
be  developed. 

The  block  diagram  in  Figure  6  shows  the  basic 
system.  The  signal  s(t)  with  unknown  amplitude  A  is 
weighted  by  the  antenna  functions  g^,  (u)  and  g^(u).  The 
channel  noise  n(t)  is  added  to  each  channel  output  to 
produce  the  observations  r^(t)  and  r2(t).  These  observa¬ 
tions  enter  the  estimator  which  produces  an  output  (an 
estimate  of  the  target  position  relative  to  boresite)  that 
is  used  as  input  to  the  angle  servo.  The  angle  servo 
then  moves  the  radar  boresite. 

The  Estimator 

A  maximum  likelihood  estimate  will  be  taken  with 
respect  to  the  observations  r^(t)  and  r»(t).  To  do  this 
a  joint  probability  density  function  of  r^(t)  and  r2(t) 
is  needed. 


Fig.  6.  Block  Diagram  of  the  Radar  Model  Showing 
the  Addition  of  the  Estimator  and  Servo 

The  observations  r^{t)  and  r2(t)  are 


r1  (t)  =  Agj,  (u )  s  (t.) 

+  n  ( t) 

(12) 

r2(t)  =  Ag^(u) s(t) 

+  n(t) 

(13) 

where  n(t)  is  additive  white  gaussian  channel  noise  and  A 
is  the  unknown  amplitude  of  the  target  return. 

For  convenience,  let  g^,  (u)  =  g^(u)  and  g^(u)  = 

g2(u) . 

The  maximum  likelihood  estimator  is  formed  by 
operating  on  the  density  function  of  the  observations. 
Since  the  observation,  r(t)  =  As ( t) g (u) +n ( t) ,  falls  under 
the  category  of  an  uncountable  infinite  number  of  points, 
it  needs  to  be  transformed  into  a  set  of  countable 
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(15) 


Since  the  noise  components,  n^,  are  independent  gaussian 
random  variables,  r^  and  r 2^  are  gaussian  and  independent. 
Their  statistics  are 


Etr^/AjU)  =  A/Eg^  ( u )  ;  Var  (r.,  ^AfU)  =  Nq/2  (16) 
E{r21/A,u}  =  A/Eg2(u);  Var(r21/A,u)  =  Nq/2  (17) 

At  this  point  only  a  single  pulse  is  being  con¬ 
sidered  so  the  second  subscript  on  the  vector  observations 
will  be  dropped. 

The  pdf  of  the  observations  is  the  joint  pdf  of 
r1  and  r2.  Since  they  are  independent,  the  joint  pdf  is 

p(rLr2/A,u)  =  p(r1/A,u)p(r2/A,u) 


P  (rir2/A,u)  ^IIN  ^ 

o 


exp  {  —  [  (r^-A/Eg^u) )  2  +  (r2~A/Eg2  (u) )  2]  } 
°  (IB) 


we  find 


To  obtain  the  maximum  likelihood  estimate  [Ref  6] 


3/3A[£n  p(r1r2/A,u)]  =  0 


3/3uUn  p(r1r2/A,u)]  =  0 
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Let  3/3ug(u)  =  g  (u) 


£n  p(r1r2/A,u)  =  ink  -  1/Nq  [  (r1"A/Eg1  (u)  )  2 
+  (r2-A/Eg2(u) ) 2] 

where  k  is  a  constant. 

Then 

3/3A[£n  p  (r1r2/A,u)  ]  =  ^-[  2  (r^^-A/Eg^^  (u)  )  (-g1  (u)  /E) 

o 

+  2(r2“A/Eg2 (u) ) (-g2 (u) /e) ]  =  0 
Solving  for  A  we  get 

~  r  g  (u)  +  r 0g0 (u) 

a  =  - v — •  <19) 

/E  (g^  (u)  +  g2  (u)  ) 

Now  to  find  u 

3/3u[ in  p  (r1r2/A,ii)  ]  =  ~[  2  (r1-Av/Eg1  (u)  ) 

o 

(-A/Egj^  (u))+2  (r2~A/Eg2  (u)  )  (-A/Eg2  (u)  )  ]  =  0 

Solving  for  a/e  yields 

•  • 
r, gn  (u)  +  r (u) 

Av^E  =  -  . - — - ; - 

gilu)gi(u*  +  52<u)g2(u) 

Substituting  equation  (.19)  for  a/e  yields 

•  • 

rlgl*U*  +  r2g2^u)  r191(«)  +  r252(u) 

g^tu)  +  g22(u)  g^ujg^u)  +  g2(u)g2(u) 
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continuing  we  get 


•  • 

^rlgi  +r2g2  (u)  (u)  +92  (u)g2  (u)  ] 

=  [r.^  (u)  +r2g2  (u)  ]  [g^2  (u)  +g22  (u)  ] 

2  *  *  * 
rlg1  (u)gl<u)  +  r^g^ (u)g2(u)g2(u)  +  r2g1 (u)g2 (u) g  (u) 

2 

+  r2g2  (u)g2<u) 

2  2  2 

=  r1g1  (ujg^u)  +  rxg2  (u)g1(u)  +  r2g1  (u)g2(u) 

+  r2g22(u)g2(u) 

•  • 

r1g2(u) [gx (u) g2(u) -g1 (u)g2 (u) ] 

•  • 

+  r2g1(u)[g1  (u)  g2 (u) -gL (u) g2 (u) ]  =  0 

A 

which  for  u  =  u  results  in 

rlg2(u)  “  r2gl(u)  =  0  (20) 

Equation  (20)  is  the  general  error  equation  for  a 
one -dimensional  Amplitude-Comparison  Monopulse  Radar.  This 
equation  will  be  used  to  steer  the  radar  boresite  to  the 
target.  To  see  this,  recall  that  g^(u)  =  cos  u  and 
g2(u)  =  sin  u.  Substituting  these  into  equation  (20) 
gives 

/S  /N 

r1  sin  u  -  r2  cos  u  =  e  (21) 
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The  expected  value  of  equation  (21)  is 

A  A 

E{r^  sin  u  -  cos  u} 

A  A 

=  A/E [cos  u  sin  u  -  sin  u  cos  u] 

A 

e  =  A/E  sin(u  -  u)  (22) 

This  represents  the  error  equation  for  a  no  noise  case. 

A. 

Clearly,  if  a  =  u  the  error  is  zero.  This  means  that  the 
radar  boresite  is  pointing  directly  at  the  target.  When 
the  WGN  is  considered,  however,  the  error  is  never  zero 
and  the  error  equation  is  random. 

/v  A 

e  =  r1  sin  u  -  r2  cos  u  (23) 

The  Discriminator  Curve 

If  equation  (22)  is  normalized  with  respect  to 
the  return  signal  amplitude  and  energy  it  becomes 

e  =  sin  (Au)  (24) 

A 

where  Au  =  u  -  u.  This  equation,  which  is  the  mean  of 
equation  (23),  describes  a  discriminator  curve.  The  dis¬ 
criminator  curve  can  be  used  to  obtain  a  value  for  Au 
which  is  a  function  of  the  target  angle  0  (see  Figure  7) . 

For  each  value  of  e  on  the  discriminator  curve 
there  is  a  unique  corresponding  value  for  Au  as  long  as  the 
linear  portion  of  the  curve  is  used. 
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Fig.  7.  Discriminator  Curve 

Equation  (23)  is  also  a  discriminator  curve.  How¬ 
ever,  since  r^  and  r 2  are  random  variables,  a  different 
curve  is  produced  for  each  value  of  and  r This  pre¬ 
vents  the  use  of  equation  (23)  to  design  a  discriminator. 
Equation  (24),  being  the  average  or  mean  of  (23),  is  more 
practical  for  a  discriminator  design  since  it  is  unchanging 
with  each  observation.  Both  equations  (23)  and  (24)  are 
used  to  produce  a  Au  to  steer  the  radar  boresite.  Equa¬ 
tion  (23)  produces  the  error  value  e  which  contains  the 
target  position  information  plun  a  noise  component.  Equa¬ 
tion  (24) ,  or  the  discriminator  curve  in  Figure  7,  uses 
the  value  of  e  from  equation  (23)  to  produce  a  Au  that 
does  not  include  the  noise  component.  When  Au  is  sent  to 
the  angle  servo  to  move  the  boresite,  a  tracking  error 
develops . 
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This  can  best  be  seen  by  an  example, 
target  position  relative  to  boresite  be  1.0°.  Ideally, 
the  error,  e,  would  produce  a  Au  on  the  first  return  pulse 
that  corresponds  to  1.0°.  However,  due  to  the  additive 
WGN,  let  Au  from  Figure  7  correspond  to  0.9°.  This  value 
is  sent  to  the  angle  servo  to  move  the  radar  boresite  and 
the  resulting  new  target  position  relative  to  boresite 
is  0.1°.  This  is  also  the  tracking  error.  On  the  next 
pulse  an  error  is  produced  tha..  corresponds  to  a  target 
angle  of  0.1°  plus  a  noise  component.  Each  time  (i.e., 
observation) ,  the  boresite  is  moved  closer  to  the  target 
until  a  steady-state  error  is  developed.  As  will  be  shown 
later,  the  steady-state  error  can  be  formed  into  an  MSE. 

To  quantify  the  threshold  of  the  radar  system  the. 
MSE  needs  to  be  lower  bounded.  The  lower  bounds  will 
reveal  the  minimum  the  MSE  can  be  for  the  system  and  also 
help  determine  the  SNR  level  at  which  breaklock  or  the 
threshold  effect  occurs.  The  two  bounds  to  be  considered 
are  the  Cramer-Rao  and  Ziv  Zakai  bounds  [Refs  6;  12]. 

The  Ziv-Zakai  (ZZ)  Bound 

The  ZZ  bound  as  developed  in  references  [8;  12] 
is  stated  as  follows: 

e2(6)  >  max  sin2 0  P  (-0,9)  (25) 

“  o<e<e  e 

- p 
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where  0  is  the  target  angle  relative  to  boresite  and  0 
is  defined  by 


P 


sin  0 

P 


sin  0  +  sin  0 

_ max _ 

2 


The  probability  of  error,  Pe(-0,0)  is  determined  from 
detection  theory  and  denotes  the  error  probability  of  the 
best  procedure  for  deciding  whether  a  target  is  at  -0  or 
0  when  it  is  known  that  a  target  is  at  one  of  these  posi¬ 
tions  and  each  has  equal  possibilities  [Refs  6;  8]. 

The  probability  of  error  can  be  described  as 

f 

P(e)  =  P{Hq}P{A>0/Hq}  +  P{H1>P{JI<0/H1}  (26) 

where  HQ  and  are  two  observation  hypotheses  and  &  is 
the  log-likelihood  ratio  as  defined  in  reference  [6] . 

For  the  minimum  probability  of  error  case,  P{HQ}  =  PtH^} 
and  therefore  the  detection  threshold  =  1  [Ref  6] . 

To  begin  with,  the  two  hypotheses  are 

Hq:  -  A/Egj^  (uQ)+n1 ;  r  2  =  A/Eg2(uQ)  +  n2  (27) 

Hl:  rl  =  A/Eg^  ^ui)+n1 ;  2  =  A/Eg2(u1)  +  n2  (28) 

Recall  that  u  =  kQx  sin  0.  The  probability  densities  are 

e(rir2/Ho)=  (Hih  e*p{ir[<ri-  A^fUo))2 

o  o 

+ ‘(r2-A/Eg2(uo))2]}  (29) 
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(30) 


P(rlr2/Hl)  =  <nN-)  (u2) )  2 

o  o 

+  (r2-A/Eg2(u1))2]} 

The  likelihood  ratio  is  formed  from  the  ratio 
of  (30)  to  (29) 

exp^Ur^A/Eg^i^))2  +  (r2»A/Eg2  (v^)  )  2] }  H1 

exp{~[  (r1-A/Eg1(uQ)  )  +  (r2~A/Eg2 (uq) )  ) }  * 

o  n2 

The  log-likelihood  ratio  is 

JlnMr)  =  rx  [g1  (ux) -g^^  (uQ)  ]  +  r  2  tg2  (uL) -g2  (uq)  ] 

H1 

+  ^~lg12(u0)-g12(u1)+g22(uo)-g22(u1))  J  0 

Ho 

(31) 

The  log -likelihood  ratio,  denoted  l,  has  statistics  as 
follows : 

EU/Hq}  »  A/E[g1(u1)g1(uo/-g12(uo)+g2(uo)g2(u1)-g22(uo)) 
A/E ,  2.  .  2,  .  .  2,  >  2 , 

+  “2“ <Vgl  <V+g2  <Uo)'92  U1 

=  A/E[g1(uo)g1(u1)+g2(uo)g2(u1)] 

-  ^(g12(u0)+g12(u1)+g22(u0)+g22(u1n 

-  ~-r  Eo(1-po'  ,32) 
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where 


:o  =  912(uo)  +  912(“1>  +  922(uc)  +  922(u1) 


po  =  2/E0lg1(uo)g1(u1)  +  92(u0l92<uin  (32a) 


Similarly, 


e{  a  /h1>  = 


EoU't>o) 


The  variance  of  i  is 


Var  U/Hq)  =  E{  (l -I)  2} 

=  N0/2[gl2(uo)+gl2(ul)+g22(uo)+g22(ul) 

~  2(g1(uQ)g1(u]L)+g2(uo)g2(u1) )  ] 

=  No/2  E0(l-p0)  (34 

From  these  statistics,  the  pdfs  of  Z  given  each 
observation  hypothesis  Hq  and  are 


pU/H.1  =  (-i-2)’>  expt-^n^  ^  E0U-p0))2]}  (35) 

u  23la.  2a  „ 


pU/H.)  =  ( — ^-5)  exp{—  ji  U-  ^  EQ(l-p0)  )  21)  (36) 

2noj  2ot 


where 


“n  "  V2  Eo,1'po) 


For  this  situation,  P{£>0/Hg}  =  P{&<0/H^}  and  the  proba¬ 
bility  of  error  then  becomes 


P(e)  =  1/2PU>0/Hq}  +  1/2  P{2.<0/H1)  =  pU>0/Hq} 


-/%«*/  h0 


)ai 


exp{-^[<Ji  +~  E  (1-P  ))2]}  (37) 

2a  ^  to  o 


This  integral  can  be  simplified  with  a  change  of  variables 
by  letting 


U+  At/E/ 2  E  ( 1- P  )  ) 
o  o 


Therefore 


Substituting  back  into  (37)  yields 

r°°  i  "x2/2 

P(e)  =/  — - —  e  dx  (38) 

Jx'  /~2n 


where 

X*  =  A/E[E  / 2N  (1-p  )  J** 
o  o  o 

Equation  (38)  can  be  recognized  as  the  Q-function  (Ref  6] . 
The  probability  of  error  is  now 

E  E 

P(e)  =  Q(x')  =  QtA(-2^£-(l-P0)  )**)  (39) 
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If  symmetry  in  the  error  is  assumed  (i.e.,  uq  =  -u^  =  u) 
then 

.2  2  2  2 
=  sin  u  +  sin  u  +  cos  u  +  cos  u  =  2 
o 

and  from  equation  (32a)  with  =  -u  =  uQ 

2  2 

PQ  =  (-sin  u  +  cos  u)  --  cos  2u 

The  variable  P  is  known  as  the  correlation  coefficient 
o 

[Ref  6]  and  its  result  is  consistent  since  it  must  vary 
between  ±1.  Equation  (39)  becomes 

P  (e )  =  Q[A(E/Nq(1  -  cos  2u))Js]  (40) 

Having  determined  the  probability  of  error,  the  ZZ 

bound  can  now  be  derived  for  the  MSE  of  the  radar  model. 

If  we  substitute  u  =  k  x  sin  6  and  let  k  x  =  1,  then  for 

o  o 

an  MSE  of  e2(0)  the  lower  bound  is  [Ref  12] 
e2(u)  ^  u2  Q [A(E/Nq(1-cos  2u))is] 

and  then 

e2 (0 )  >  sin2  0  Q [ A ( E/Nq ( 1-cos ( 2  sin  0)))^]  (41) 

Two  properties  can  be  seen  from  equation  (41) . 

First,  for  large  SNR  (Eq/No)  the  argument  of  the 

Q-function  will  be  large,  giving  a  small  value  for  Q(x) . 

This  means  a  small  probability  of  error  and  hence  a 

smaller  MSE.  This  is  intuitively  consistent  since  we  would 

2 

expect  to  track  well  at  high  SNR.  As  0  increases,  f  (Q) 


25 


increases,  which  tells  us  that  the  M5E  becomes  larger  when 
we  are  not  pointing  at  the  target.  So  in  order  to  keep 
the  MSE  small  the  boresite  must  be  kept  on  the  target. 

The  Cramer -Rao  Bound 

In  reference  [2],  it  was  shown  that  for  a  one¬ 
dimensional  Amplitude-Comparison  Monopulse  Radar,  the  CR 
bound  for  the  variance  of  the  estimated  target  position, 

A 

Var  (e  ) ,  was 


Var (0  )  >  5 -  (42) 

A  N  (gl  (u)+g2  (u)) 

’o 

•  2  2  i  *2 

For  this  case  g^  (u)  =  sin  u  cos“  0  and  g2  (u)  = 

2  2 

cos  u  cos  0  so  (41)  becomes 

Var  (0  )  >  - ?  pp1 - 5— 7r=-p - - - a  (4  2a) 

A  ~  cos  0  A  SNR  cos  0 
No 

*  2  *  2 

The  factor  (g^  (u)  +  g2  (u) )  can  be  shown  to  be 
the  inverse  of  the  beam-width  [Ref  2] . 

The  CR  bound  states  that  a  narrow  beam-width  will 
ensure  a  small  variance.  But  the  CR  bound  does  not  account 
for  the  fact  that  when  the  threshold  effect  and  breaklock 
occur,  the  target  is  outside  of  the  radar  beam-width  and 
the  variance  of  the  error  is  very  large.  At  this  point 
(i.e.,  very  low  SNR)  a  narrow  beam-width  no  longer  guar¬ 
antees  a  small  error  variance.  The  bound  becomes 
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unreliable  and  cannot  accurately  lower  bound  the  MSE  at 
low  SNR  levels. 

The  Tracking  Loop 

To  gain  insight  into  where  the  breaklock  point 
occurs  it  is  necessary  to  design  a  tracking  loop  which 
uses  equation  (22)  to  steer  the  antenna  boresite  toward 
the  target.  Figure  P  shows  a  block  diagram  of  the  mono 
pulse  system  which  will  be  simulated. 


Servo 


Fig.  8.  Block  Diagram  of  AC  Monopulse  Radar  Showing 
the  Estimator  being  used  to  Track  a  Target 

In  Figure  8,  s^ ,  s 2,  n^  and  n^  are  the  representa¬ 
tions  of  the  signal  and  noise  components  respectively, 
where  s^  =  A/Eg^fu)  and  s2  =  A^Eg^u)  .  The  observation 
r^  =  s^ +  n^  and  r 2  =  s2  +  n2  are  random  numbers.  For 
simulation  purposes  they  are  determined  by  using  a  random 
number  generator  to  generate  n^  and  n2  which  are  then 
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added  to  and  s 2  respectively.  The  antenna  functions 
are  determined  by  noting  that  equations  (10)  and  (11)  can 
be  written  as 


g^(u)  =  cos  (sin 

9) 

(10) 

g2(u)  =  sin  (sin 

0) 

(11) 

where  k  x  =1  (this  means  that  the  distance  between  the 
o  o 

two  point  sources  is  small  and  has  been  chosen  so  that 

k  x  =  1) . 
o  o 

The  target  position  relative  to  the  boresite  is  0. 
Since  this  is  a  computer  simulation,  the  actual  target 
position  is  inserted  into  the  antenna  functions  to  provide 
the  amplitude  weighting.  With  this  information,  r^  and 
r2  are  calculated  and  fed  into  the  estimator.  The  esti¬ 
mator  to  be  used  is  equation  (23)  which,  restated,  is 

A  A. 

e  =  r^  sin  (sin  6)  -  r2  cos  (sin  6)  (23) 

The  error,  e,  is  calculated  using  the  estimated 

A 

target  position  8,  and  the  random  numbers  r^  and  r 2  (which 
are  calculated  using  the  actual  target  position  8) .  The 
target  position,  estimated  and  actual,  is  always  relative 
to  boresite.  If  we  assume  that  the  radar  boresite  is 

A 

pointing  to  the  estimated  target  position  8,  then  for  each 

A 

observation  8=0  and  the  resulting  error,  e,  contains 
information  pertaining  to  the  actual  target  position. 
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The  error  e ,  which  now  contains  the  target  posi¬ 
tion  0  and  a  noise  component  is  fed  to  the  discriminator. 
Recall  that  the  mean  of  the  error  e  is  equation  (22) 
rewritten  here  as 

_  A 

e  =  sin  (u  -  u) 

A 

=  sin  (sin  6  -  sin  0)  (43) 

If  a  3  dB  beam-width  of  3  degrees  is  assumed, 

A0  ,  or  the  difference  between  boresite  and  the  target 
will  be  a  maximum  of  ±  1.5  degrees.  For  small  values  of  0, 
sin  0^0,  therefore 

e  =  sin  (0-0)  (44) 

If  we  assume  that  the  boresight  is  pointing  at  the  esti¬ 
mated  target  position,  then 

e  =  -  sin  e  (45) 

Equation  (45)  can  be  further  reduced  to 

e  =  -0  (46) 

For  the  range  of  0  chosen,  the  actual  discriminator  curve 
is  shown  in  Figure  9. 

Because  of  the  assumptions  made,  0  in  equation 
(46)  is  the  same  as  A0 ,  the  distance  in  degrees  the  bore- 
site  must  move  to  be  pointing  at  the  target. 
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Fig.  9.  Discriminator  Curve  Showing  Maximum 
Values  of  e  and  6  or  A6 

As  stated  before,  the  error  is  transformed  into  a 
corresponding  A0  by  use  of  the  discriminator  curve.  The 
error  e,  fed  to  the  discriminator  is  random  due  to  the 
additive  channel  noise.  The  discriminator  curve  used 
represents  the  error  without  noise;  therefore,  the  corres¬ 
ponding  A0  produced  is  also  random. 

The  next  step  in  the  tracking  loop  is  to  input 
the  target  position  0,  or  A9 ,  into  the  angle  servo  to 
steer  the  boresite  to  the  target.  The  value  used  for  A0 
will  be  constant  during  each  observation  interval;  there¬ 
fore,  it  can  be  modeled  as  a  step  input  to  the  servo.  The 

i 

amplitude  of  the  input  to  the  servo  after  the  next  observa¬ 
tion  will  reflect  the  change  in  the  radar  boresite  as  it 
moves  closer  to  or  away  from  the  target.  The  input  becomes 
a  series  of  step  functions  whose  amplitudes  are  proportional 
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to  the  distance  the  boresite  must  move  (in  either  direc¬ 


tion,  ±0)  to  the  target. 

The  output  of  the  servo  will  be  delayed  due  to 
the  inertia  of  the  antenna  and  the  drive  motor,  but  it 
needs  to  be  approximately  equal  to  its  input.  The  servo 
used  should  therefore  provide  a  zero  steady-state  error 
for  a  step  input. 


Fig.  10.  Angle  Servo  Feedback  Control  Loop 

Assume  the  angle  servo  has  a  Jeedback  control  loop 
as  shown  in  Figure  10.  The  transfer  function  is 


K 

A0out  ^ 

s(s+w) 

A0 

A0 


out 

in 


2 

s 


K 

+  sw  +  K 


Equation  (47)  can  be  factored  into 


G  ( s) 


K 

(s+a) (s+g ) 


(47) 


(48) 
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Using  (48)  as  the  servo  transfer  function  and  assuming 
the  input  to  the  servo  to  be  a  step  function  A0(s)  =  AG/s, 
the  servo  output  becomes 


A8/s  •  G ( s)  = 


A6K 


s ( s+a )  ( s+$  ) 


which  has  an  inverse  Laplace  transform  of 


40K  A9K(Be~at  -  ae~8t) 
a3  a3  (a -3 ) 


(49) 


(50) 


If  we  let  a  =  -a  +  joj ,  3  =  ~a  -  jw  and  K  =  a  3»  then  after 
some  manupulation,  equation  (50)  becomes 

A0 ' (t)  =  A8 [1-e  at(—  sin  wt  -  cos  ut) ]  (51) 

U) 

The  values  of  a  and  w  are  obtained  from  the  poles  of  the 
transfer  function.  The  pole  positions  are  chosen  so  that 
a  value  of  0  will  damp  the  servo  system  response  such  that 
it  will  settle  out  without  four  system  time  constants. 

This  settling  time  will  give  an  output  within  5  percent  of 
the  input  value.  The  settling  time  is  chosen  to  be  a  multi¬ 
tude  of  the  PRI  of  the  radar.  The  PRI  used  for  the  computer 
simulation  is  1/200  pps  or  5  ms. 

Once  an  output  A0 '  is  obtained  from  the  servo, 
it  is  used  to  steer  the  radar  boresite.  This  is  simulated 
by  changing  the  antenna  functions  g^(u)  and  g2(u)  tc  ref^ect 
the  new  target  position.  If  there  were  no  noise,  the 
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boresite  would  point  directly  at  the  target  after  one 
pulse.  Because  of  channel  noise,  however,  the  boresite 
will  point  close  to  the  target  after  each  pulse  but  never 
directly  on  it.  The  new  (or  adjusted)  antenna  functions 
produce  a  new  error  e  which  in  turn  produces  a  new  A0 . 

This  loop  continues  as  the  radar  boresite  tracks  the  target 
within  some  small  error  (MSE) . 

To  begin  the  problem,  a  stationary  target  is 
tracked  and  the  MSE  is  calculated.  The  stationary  target 
is  a  simple  case  since  in  the  routine  the  target  position 
is  initialized  within  the  3  dB  beam-width  and  tracking  it 
consists  of  maintaining  the  boresite  on  it. 

To  consider  a  moving  target,  a  relationship  between 
the  target  position,  velocity,  and  range  is  needed. 

Figure  11  shows  the  geometry  of  the  situation.  The  dis¬ 
tance  .the  target  moves  in  the  positive  x  direction  is 
V^t ,  where  VT  is  the  target  velocity.  Rq  is  the  range  to 
the  target  at  the  point  of  acquisition  and  is  the 
changing  range  to  the  target.  Note  that  the  target  range 
RT  and  the  SNR  are  inversely  proportional.  The  changing 
SNR  can  be  accounted  for  by  using  the  radar  range  equa¬ 
tion  [Ref  1] . 


.  P  G  2  X2  o 

SNR  R4  =  -  R  — ^ - - 

<4II)  N 


(52) 
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Fig.  11.  Geometry  Showing  the  Relationship  Between 
Target  Position,  Velocity  and  Range 


where 


=  range  to  the  target 
=  power  received 
=  receiver  gain 
=  wavelength  of  carrier 
=  target  cross-section 
=  internal  receiver  noise  power 


For  the  problem  at  hand,  the  right  side  of  equa 
tion  (52)  remains  constant.  If  we  let  SNRq  and  Rq  be 
the  initial  SNR  and  target  range  and  SNRT  and  R^  be  the 
new  SNR  and  target  range  we  have 


SNR  R 
o  o 


PR  S2  ^  °T 

( 4  n ) 3  n 


snrt  Rt 


PR  Gr;  *2  °T 

(4n> 3  n 
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Therefore 


SNRt  ^  =  SNR0  R04 


and 


snrt 


SNR 

o 


R  4 


(55) 


Equation  (55)  can  account  for  the  changing  SNR  as  the 
target  range  R^  increases  for  a  moving  target. 


tan  eT 


eT(t)  =  tan_1(VTt/Ro) 


(56) 


With  a  relationship  for  0T  as  a  function  of  velocity, 
time,  and  range,  the  target  position  can  be  accounted  for 
in  the  computer  simulation. 

The  Error 

The  tracking  error  for  the  stationary  and  non¬ 
stationary  case  is  the  difference  between  the  radar  bore- 
site  and  target  position. 

ET  =  eB  -  °T  (57> 

The  simulation  will  take  50  observations  which  could  con¬ 
sist  of  every  pulse  or  every  several  pulses  depending  on 
the  settling  time  chosen  for  the  angle  servo.  Each  set 
of  50  observations  is  considered  a  run.  In  each  run  a 


"seed"  is  used  to  generate  a  vector  of  random  numbers 
which  are  used  to  determine  n^  and  n2  which  in  turn  are 
used  to  calculate  r.^  and  r^.  A  different  seed  is  used  for 
each  separate  run  and  a  new  set  of  random  data  is  generated 
If  equation  (57)  is  squared  for  each  observation  and  then 
averaged  over  the  number  of  runs,  an  approximation  of  the 
MSE  is  obtained.  For  an  example  see  Figure  12.  In 
Figure  12a  the  value  at  observation  4  is  summed  over  all 
the  runs  and  then  divided  by  the  number  of  runs  to  obtain 
the  MSE  for  that  observation  which  is  plotted  in  Figure  12b 
If  this  is  done  for  each  observation,  the  approximate  MSE 
is  obtained  for  a  single  SNR  level.  Recall  that  an  obser¬ 
vation  can  be  one  pulse  or  a  multiple  of  pulses. 

The  threshold  or  point  of  breaklock  can  be  deter¬ 
mined  by  observing  when  the  MSE  increases  continually  per 
observation.  As  the  SNR  is  decreased,  the  plots  of  the 
MSE  will  look  similar  to  Figure  13. 

Since  the  beam-width  is  ±  1.5  degree,  an  error  e 
producing  a  A0  greater  or  less  than  ±  1.5  degrees  will  be 
automatically  limited  to  +1.5  degrees.  The  servo  will  then 
move  the  boresite  its  maximum  limit  and  another  observa¬ 
tion  will  be  made.  If  the  error  developed  at  this  obser¬ 
vation  produces  a  A9  again  greater  than  the  maximum, 

A0  will  automatically  be  limited  again  to  ±1.5  degrees. 

The  servo  will  not  be  able  to  keep  up  with  the  moving 
target  to  be  able  to  track  it.  This  is  breaklock  and  can 
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Pulses 


b 

Fig.  12.  (a)  .ample  Plots  of  the  Error  Squared  per 

Observation  for  One  SNR  Level;  (b)  Example  Plot  of 
the  MSE  per  Observation  for  One  SNR  Level 


MSE 


1 _ _ _ ► 

Pulses 

Fig.  13.  Sketch  of  MSE  Becoming  Unbounded 
at  Breaklock 

be  caused  by  a  number  of  factors.  The  target  velocity 
can  be  fast  enough  at  a  certain  range  to  cause  breaklock, 
the  servo  could  be  too  slow  to  "catch  up"  to  the  target, 
or  the  SNR  can  degrade  to  a  point  where  the  threshold 
effect  occurs.  The  velocity,  range,  and  servo  will  be 
adjusted  to  prevent  breaklock  except  at  a  low  SNR. 
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III.  The  Computer  Program 


The  programs  for  the  stationary  and  nonstationary 
target  are  straightforward.  The  only  subroutines  used  were 
two  IMSL  subroutines,  GGNML  and  MDNOR. 

The  program  was  designed  so  that  the  radar  param¬ 
eters  such  as  PRI,  servo  settling  time,  and  boresite  posi¬ 
tion  can  be  initiated  at  the  beginning.  The  initial  target 
position  must  be  within  ±1.5  degrees  of  the  boresite  posi¬ 
tion  . 

The  IMSL  subroutine  GGNML  is  called  to  generate  a 
vector  of  random  numbers  to  represent  the  noise  components 
n^j^  and  n2  of  the  observations  r^  and  r2 .  A  vector  of  100 
random  numbers  is  generated  and  two  numbers  arc  used  for 
each  observation.  The  program  is  set  up  for  50  observa¬ 
tions  but  can  be  adjusted  for  any  amount  provided  the  random 
number  vector  is  also  adjusted.  The  random  numbers  r^  and 
r2  are  then  calculated  using  equations  (14)  and  (15) .  The 
value  for  6  is  calculated  by  subtracting  the  boresite  posi¬ 
tion  from  the  target  position.  It  is  then  updated  as  the 
boresite  is  moved.  Once  an  error  is  calculated  from 
equation  (23),  it  is  converted  to  a  A0  by  use  of  equation 
(24).  This  A0  is  used  as  input  to  the  servo  equation  (51). 
The  output  of  equation  (51)  is  used  to  move  the  antenna. 
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This  is  simulated  by  changing  g  (u)  and  g^ (u)  to  reflect 
the  new  difference  between  radar  boresite  and  the  target. 

For  a  moving  target,  equation  (54)  is  used  at  this  point 
to  calculate  the  new  target  position.  The  SNR  will  change 
as  the  range  from  the  radar  to  the  target  increases.  This 
is  accounted  for  in  the  program  for  the  moving  target.  The 
boresite  position  is  updated  by  adding  the  servo  output  to 
the  old  boresite  position. 

All  of  the  above  are  accomplished  for  one  observa¬ 
tion.  After  this  is  done  the  tracking  error  is  calculated 
by  subtracting  the  target  position  from  the  boresite  posi¬ 
tion.  This  value  is  squared  and  stored  in  a  50  column  by 
15  row  array.  The  ZZ  and  CR  bounds  are  also  calculated  at 
this  point  and  stored  in  a  similar  array.  After  15  runs, 
the  columns  of  the  array  are  summed  and  divided  by  the 
number  of  rows  (15)  to  form  the  mean- square  error  per 
observation.  The  MSE  is  then  plotted  for  50  observations. 
The  entire  process  is  then  repeated  for  a  lower  SNR  level. 

The  SNR  level  was  varied  from  35  dB  to  -10  dB. 

The  target  velocity  used  was  609.6  meters/ second.  The 
servo  settling  time  used  was  100  milliseconds.  The  acquisi¬ 
tion  range  was  6000  meters. 

To  calculate  the  ZZ  bound,  the  IMSL  subroutine 
MDNOR  was  called  to  evaluate  the  Q-function  used  to  deter¬ 
mine  the  probability  of  error. 
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Fig.  14.  Flow  Diagram  for  Stationary/Nonstationary 

Target  Tracking 
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Fig.  14 — Continued 


IV.  Results  and  Conclusions 


General 

The  computer  program  was  run  for  the  stationary 
and  the  nonstationary  target  with  the  SNR  varied  from 
35  dB  to  -10  db.  For  each  SNR  level,  15  runs  were  made 
taking  50  observations  each.  The  MSE  per  observation  was 
computed  and  plotted.  In  addition  to  the  MSE,  the  Cramer- 
Rao  and  Ziv-Zakai  bounds  were  calculated  per  observation 
and  plotted  (see  Appendix  B and  C) . 

The  object  was  to  determine  the  point  of  breaklock 
and  how  the  lower  bounds  can  be  used  to  quantify  it. 

Stationary  Target 

As  the  SNR  level  was  decreased,  the  MSE  increased. 
It  became  increasingly  difficult  for  the  boresite  to  stay 
on  target.  The  MSE  became  unbounded  between  -7  dB  and 
-10  dB.  Breaklock  is  assumed  to  have  occurred  here. 

Cramer-Rao  Bound 

As  Figure  15  shows,  the  MSE  at  7  dB  is  greater 
than  the  constant  level  of  the  CR  bound.  As  the  SNR  was 
decreased  the  CR  bound  became  very  inaccurate.  In 
Figure  16  the  MSE  is  below  the  CR  bound  in  some  places. 

It  appears  that  the  CR  bound  is  the  average  of  the  actual 


43 


0NRRY  TARGET 


Stationary 


STRTIONRRY  TARGET  (-7  DB  ) 


00*0fr 


00  *0E 


00*03 


00*01 


QNnoa  oyy-y3wuy: 


00*0 


Stationary  Target 


MSE.  But  according  to  the  definition  of  the  lower  bound, 
the  MSE  must  always  be  greater  than  or  equal  to  the  lower 
bound.  Fiqure  17  shows  the  MSE  compared  to  the  CR  bound 
at  -20  dB .  The  MSE  is  completely  below  the  CR  bound. 

The  Ziv-Zakai  Bound 

The  ZZ  bound  for  the  stationary  target  tends  to 
follow  the  general  trend  of  the  MSE.  Note  the  large  magni¬ 
tude  difference  between  the  ZZ  bound  and  the  MSE  (see 

2 

Figure  18) .  The  ZZ  bound  is  the  product  of  the  sin  0  and 
the  P(e).  Neither  is  ever  greater  than  1;  therefore,  the 
ZZ  bound  is  never  greater  than  1.  The  importance  of  the 
ZZ  bound  in  this  case  therefore  is  that  it  shows  the  trend 
of  the  MSE.  Note  in  Figure  18  that  the  ZZ  bound  at  many 
points  has  the  same  peaks  and  valleys  as  the  MSE.  In 
Figure  19  the  MSE  becomes  unbounded.  This  trend  is  also 
shown  in  the  ZZ  bound. 

Moving  Target — Cramer-Rao  Bound 

Figures  20  and  21  show  the  MSE  as  compared  to  the 
CR  bound.  Breaklock  occurs  below  5  dB  when  the  MSE  becomes 
unbounded.  The  CR  bound  does  not  show  any  trend  toward 
breaklock  even  at  a  very  low  SNR  (Figure  21) .  The  CR  bound 
is  generally  constant  with  a  slight  upward  trend  toward  the 
end  of  the  run.  This  is  due  to  the  SNR  decreasing  as  the 
target  moves  horizontally  past  the  radar.  The  upward 
movement  of  the  bound  is  not  related  to  the  actual  MSE. 
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Fig.  18.  Stationary  Target  (1  dB) 
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Fig.  20.  Moving  Target  (5  dB) 
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As  the  SNR  is  decreased,  the  CR  bound  tells  us  that  the 
MSE  is  increasing  but  gives  no  indication  when  breaklock 
occurs . 

Ziv-Zakai  Bound 

As  in  the  stationary  case,  the  ZZ  bound  follows 
the  trend  of  the  actual  MSE.  The  ZZ  bound  becomes  unbounded 
as  the  MSE  becomes  unbounded  (Figures  22  and  23). 

Among  other  reasons,  the  breaklock  point  occurs 
because  of  SNR  level,  the  target  range  and  speed,  and  the 
angle  servo  response  time.  The  exact  level  of  SNR  for 
breaklock  is  not  important  here  since  it  would  depend  on 
the  parameters  of  an  actual  (not  simulated)  random  system. 
What  is  important  is  that  if  these  parameters  were  taken 
into  account,  the  breaklock  or  threshold  of  an  actual  sys¬ 
tem  could  be  predicted  by  the  ZZ  bound. 

Ziv-Zakai  Bound  Compared  to 
Cramer-Rao  Bound 

It  is  interesting  to  observe  why  the  ZZ  bound 
follows  the  actual  MSE  closer  than  the  CR  bound.  The  CR 
bound  as  stated  in  Chapter  II  is 

Var (0)  >  - - — 5 —  (4  2) 

SNR  cos  0 

The  angle  0  is  the  difference  between  the  bore- 
site  and  target  position  on  the  radar  beam  (i.e.,  the 
error).  It  is  limited  to  ±  1.5  degrees  since  the  3  dB 
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beam-width  is  3  degrees.  The  cosine  of  a  very  small  angle 
is  approximately  1  so  the  CR  bound  is  essentially  con¬ 
stant  for  each  SNR  level  (i.e.,  it  is  not  influenced 
heavily  by  the  error)  .  On  the  other  hand,  the  ZZ  bound  i<= 

e  2  (9  )  >_  sin2  0  Q  [A (E/Nq (1-cos  2U))*5]  (41) 

where  u  =  sin  0  and  6  is  the  target  position  relative  to 
boresite . 

The  probability  of  error  is  related  to  the  SNR 

level  and  is  small  for  lage  SNR  and  approaches  1  for  small 

SNR.  The  factor  that  allows  the  ZZ  bound  to  follow  the 

2  2 

general  trend  of  the  actual  MSa  is  sin  6  or  G  . 

Since  0  is  the  error  between  the  boresite  and 

.  2  2 
target  and  is  limited  to  +  1.5  degrees,  sin  0  =0  which 

can  be  related  to  the  mean-square  error.  Hence  the  ZZ 

bound  is  more  prone  to  follow  the  actual  MSE.  Because  of 

the  tremendous  magnitude  difference  between  the  ZZ  bound 

and  MSE,  the  ZZ  bound  does  not  give  a  good  indication  for 

a  lower  bound  on  the  MSE.  It  could  be  possible,  however, 

to  scale  the  ZZ  bound  to  the  actual  MSE  and  not  only  be 

able  to  predict  its  trend  but  a  more  accurate  and  realistic 

lower  bound.  More  research  would  be  required  to  develop  a 

scaling  factor. 

The  reason  that  9  is  limited  to  ±1.5  degrees  is 
because  it  has  been  assumed  that  there  are  no  side  lobes 
to  the  radar  or  that  they  are  being  neglected.  Once  the 
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target  leaves  the  3  dB  beam-width  it  can  no  longer  be 
tracked..  However,  if  there  exists  a  significant  level  of 
sidelobes,  then  6  could  be  limited  to  the  angle  of  these 
sidelobes. 

As  a  final  conclusion  it  is  necessary  to  observe 
what  information  the  two  bounds  give  as  to  how  to  keep  the 
MSE  small.  In  reference  [2]  it  was  concluded  that  the  CR 
bound  states  that  the  variance  is  small  for  large  SNR  and 
narrow  beam-width.  It  is  smallest  when  the  boresite  is 

kept  on  the  target.  This  can  be  seen  from  equation  (42) . 

2 

The  cos  0  term  is  1  when  0=0  which  means  the- 

.  2 
boresite  is  on  the  target.  The  cos  0  decreases  as  0  becomes 

larger,  therefore  increasing  the  variance.  The  ZZ  bound 

has  the  same  characteristics.  The  smaller  0  is  the  smaller 

the  MSE.  The  contribution  of  the  probability  of  error  is 

mainly  due  to  the  SNR  level. 

In  communication  theory,  the  correlation  coefficient, 
p.,  is  used  to  measure  the  similarities  between  two  signals. 

If  both  signals  are  exactly  the  same,  p  =  1  and  if  they  are 
opposite,  p  =  -1  [Ref  11] .  In  this  thesis  the  correlation 
coefficient  is  used  in  a  slightly  different  manner.  The 
correlation  coefficient  was  determined  to  be  p  =  cos  2  u 
from  equation  (32a) .  The  variable  u  is  equal  to  sin  0  and 
for  a  small  0,  sin  0~0  ;  therefore,  p  =  cos  20  .  The  angle 
0  is  the  target  angle  relative  to  the  boresite.  When  the 
boresite  is  directly  on  the  target,  0=0  and  p  =  1.  As  0 
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increases  the  correlation  coefficient  decreases.  This 
means  that  for  two  targets  symmetrically  located  on  the 
radar  beam,  the  resolution  ability  increases  respectively 
with  the  separation  angle  20 .  For  the  problem  presented 
in  this  thesis  there  is  only  one  target;  however,  there 
exists  the  possibility  of  a  false  target  at  its  symmetrical 
opposite.  By  using  the  correlation  coefficient,  the  error 
associated  with  detecting  the  target  or  its  symmetrical 
opposite  is  accounted  for  in  the  error  probability  in 
equation  (41) . 

Rec  ommenda t i on  s 

The  application  of  the  ZZ  bound  to  the  mean-square 
tracking  error  has  proven  fruitful.  It  is  well  worth  recom¬ 
mending  an  investigation  of  the  two-dimensional  case.  Once 
the  probability  of  error  is  determined,  and  a  suitable  two- 
dimensional  inequality  of  the  tchebycheff  type  is  found, 
the  ZZ  bound  _or  the  covariance  of  the  tracking  error  can 
be  determined. 

It  is  also  recommended  that  an  additional  noise 
source  be  considered  in  the  one-  or  two-dimensional  case. 

The  noise  source  can  be  a  repeater  jammer  or  chaff  cloud. 

As  the  radar  is  tracking  the  target  the  noise  source  will 
expedite  breaklock.  The  question  of  what  noise  power  is 
needed  to  cause  breaklock  and  how  can  the  ZZ  bound  help  in 
quantifying  the  MSE  should  be  addressed. 


Further  recommendations  are  to  consider  actual 
servo  parameters  for  the  radar  system  and  incorporate 
random  target  motion  to  simulate  evasive  action.  Also  an 
investigation  of  sidelobe  structure  for  a  typical  mono¬ 
pulse  radar  would  be  an  interesting  addition  to  the  track¬ 
ing  error. 
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The  Two-Dimensional  Monopulse  Radar  Model  and  the 
Two-Dimensional  Ziv-Zakai  Bound 


Introduction 

An  attempt  was  made  to  model  the  two-dimensional 
(2-D)  Amplitude-Comparison  monopulse  radar  and  derive  a 
2-D  ZZ  bound.  The  ZZ  bound  used  in  this  thesis  was  derived 
for  a  one -dimensional  (azimuth)  situation.  To  expand  che 
bound  to  2-D  is  not  a  simple  task  and,  as  will  be  shown, 
its  complexity  increases  four-fold.  The  2-D  radar  model 
is  slightly  more  complex  but  not  difficult. 

Since  it  is  simpler,  the  2-D  radar  model  is  devel¬ 
oped  first  with  its  resulting  antenna  functions.  The 
antenna  functions  derived  could  be  the  source  of  some  of 
the  difficulties  encountered  in  developing  this  approach. 
For  example,  they  are  no  longer  independent  of  each  other. 
Each  beam  (sum  and  difference)  in  each  dimension,  contains 
both  azimuth  and  elevation  information.  This  result  seems 
consistent  with  what  would  be  expected  from  a  two- 
dimensional  system.  The  problem  could  arise  from  using 
the  Gram-Schmidt  orthogonalization  procedure  to  describe 
the  observations. 

After  the  2-D  model  is  developed,  the  2-D  ZZ  bound 


is  discussed.  The  attempt  to  derive  this  bound  was 


abandoned  due  to  two  main  obstacles  encountered.  Recom¬ 
mendations  are  also  included  as  to  what  directions  may  be 
taken  to  overcome  the  obstacles. 

The  Two-Dimensional  Radar  Model 

A  2-D  Amplitude-Comparison  monopulse  radar  con¬ 
sists  of  four  antennas  placed  symmetrically  in  the  x-z 
plane.  As  in  the  one-dimensional  case,  the  antennas  are 
considered  point  sources  in  an  array.  Figure  A-l  shows  a 
simple  model  of  the  2-D  case  where  antennas  A  and  B  are 
split  into  a  sum  and  difference  and  likewise  antennas  C 
and  D.  This  configuration  produces  one  sum  beam  and  two 
difference  beams. 


The  azimuth  difference  beam  is  obtained  by  sub¬ 
tracting  the  difference  of  the  four  antennas  and  the  ele¬ 
vation  difference  beam  is  obtained  by  subtracting  their 

sums.  The  difference  beams  are 
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(A+D)  -  (B+C) 

rr 


(A-l) 


(A+B)  -  (C+D) 

n 


(A-2) 


The  sum  beam  is  obtained  by  adding  the  antenna  sums. 


(A+B)  +  (C+D) 

/T 


(A-3) 


The  sum  and  difference  beams  can  be  used  to  derive 
the  antenna  functions.  Assuming  that  the  monopulse  can  be 
modeled  as  a  2-D  array  with  four  point  sources,  the 

t 

approaching  plane  waves  can  be  described  as  (see  Figure  A-2) 


ea  = 

jKA'r 

E  e  A 
o 

(A-4) 

eb  = 

E  e^'? 
o 

(A-5) 

Ec  = 

jKc*r 

E  e  ^ 
o 

(A-6 ) 

ed  - 

jKD-r 

E  e 
o 

(A-7) 

where  E  is  a  constant, 
o 


Fig.  A-2.  Coordinate  System  Showing  Four 
Point  Sources  in  an  Array 

and 

A  A  A  /S 

r  a  x  sin  6  cos  <t>  +  y  sin  0  sin  <J>  +  z  cos  6 

A 

Performing  the  dot  product  K*r  in  equations  (A-4)  through 
(A-7)  gives 


Ea  =  E  exp[jk  (x  sin  0  cos  <f>  +  z  cos  0)]  (A 

A  O  O 

E_  =  E  exp[jk  (-x  sin  0  cos  <f>  +  z  cos  0)]  (A 

15  O  O 

E^,  =  E  exp[-jk  (x  sin  0  cos  <j>  +  z  cos  0)]  (A 

v  o  o 

Ed  =  Eq  exp[jkQ(x  sin  0  cos  <f>  -  z  cos  0)]  (A 


If  we  let  u  =  k  x  sin  0  cos  <j>  and  v  =  k  z  cos  0 
o  o 

and  use  Euler's  equation,  the  antenna  functions  can  be 


obtained  by  substituting  equations  (A-8 )  through  (A-ll) 
into  (A-l )  through  (A-3) .  If  the  magnitude  of  the  result 
is  taken,  the  antenna  functions  are  [Ref  8] 


el  =  2/5 

Eo 

cos 

u 

sin 

V 

(A-12 ) 

as  '  2/5 

Eo 

sin 

u 

cos 

V 

(A-13 ) 

l  -  2/2 

Eo 

cos 

u 

cos 

V 

(A-14 ) 

By  normalizing  these  functions,  the  final  form  of  the 
antenna  functions  becomes 


gx (u,v) 

=  sin 

u 

cos 

V 

(A-l 5) 

g2 (u,v) 

=  cos 

u 

sin 

V 

(A-16) 

g3 (u,v) 

=  cos 

u 

cos 

V 

(A-17) 

A  block  diagram  of  the  model  can  now  be  shown. 
Figure  A-3  shows  the  output  of  the  four  antennas  (target 
return)  being  weighted  by  the  three  antenna  functions. 
White  gaussian  channel  noise  is  added  to  the  weighted 
target  returns  to  form  the  observations  r^(t),  r2(t),  and 
r^(t).  The  same  assumptions  apply  to  the  WGN  as  did  in 
the  one-dimensional  case. 

In  order  to  use  the  observations  to  obtain  a 
maximum  likelihood  estimate  of  the  azimuth  and  elevation 
angles,  the  Gram-Schmidt  orthogonalization  procedure  is 
used  as  before  with  the  one-dimensional  case.  The  vector 
observations  are  therefore 
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Fig.  A-3 .  Block  Diagram  of  the  Two- 
Dimensional  Radar  Model 


rl  =  A/E  g1  (u, v)  + 
r2  =  a/e  g2(u,v)  +  n2 
r,  =  A/E  g~  (u,  v)  +  n. 


(A-18 ) 


(A-19) 


(A-20 ) 


The  mean  and  variance  of  these  vectors  are 


E{r1/A,u,v}  =  A/E  g1(u,v);  Var (r1/A,u, v)  =  Nq/2  (A-21) 
E{r2/A,u,v}  =  A/E  g2(u,v);  Var (r2/A,u, v)  =  Nq/2  (A-22) 
E{r3/A,u,v}  =  A/E  g3(u,v);  Var (r3/A,u, v)  =  NQ/2  (A-23) 

The  joint  pdf  is  again  the  product  of  the  individual  pdfs. 


p(rlr2/A/u, v)  -  (2n  N  Yj) 


exp^2lN~72)t(rl'‘A>^gl^2'f  <r2~A|/^g2)  2+(r3-A/Eg3)  2]  } 


(A-24 ) 


«*•***• 


The  maximum  likelihood  estimate  is  obtained  by  finding 
[Ref  6] 

^ip(r1r2r3/A,u,v)  =  0 
^|^p(r1r2r3/A/u, v)  =  0 

^^•p(r1r2r3/A,u,  v)  =  0 

3g.  .  3g. 

Let  -r-^tUjV)  =  g.  =  -7~(ufv) 

3u  3v 

for  notational  purposes. 

The  steps  to  obtain  the  maximum  likelihood  esti¬ 
mate  are  the  same  as  in  the  one-dimensional  case.  The 
result  for  the  A  estimate  is 


rl9l  +  r2g2  +  r3g3 

A  ,-.2  2  2 . 

/E  (g^  +  g2  +  g3  ) 

The  general  error  equation  is 

•  •  2  *  2  * 

rl(glg2g2  +  glg3g3  "  g2  gl  "  g3  gl> 

•  •  2  •  2  • 

+  r2(glg2gl  +  g2g3g3  “  gl  g2  "  g 3  g2> 

•  •  j  •  2  • 

+  r3<glg3gl  +  <32<3392  “  93  -  92  93>  =  0 


(A-25) 


(A- 26 ) 


The  partial  derivatives  of  the  antenna  functions  are 
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A  2 

£2  sin  v  cos  u(r2  sin  u  +  r^  cos  u) 

A 

-  r^  cos  u  cos  v  (A-31) 

e1  and  e2  represent  the  error  equations  for  azimuth 
and  elevation.  Note  that  both  equations  contain  azimuth 
and  elevation  information. 

The  expected  values  of  and  e2  are 

—  —  2  A 

=  A/E  cos  v  sin(u  -  u)  (A-32) 

e2  =  a/e  cos  u  sin(v  -  v)  (A-33) 

A 

With  these  results  it  can  be  seen  that  if  u  =  u  then 

/\ 

=  0  and  if  v  =  v,  £2  =  0.  If  u  represented  an  azimuth 
bearing  and  v  represented  an  elevation  bearing,  these  two 
equations  could  be  used  separately  to  steer  the  radar  in 
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each  dimension.  However,  recall  that  u  =  kQ  x  sin  8  cos  <j> 
and  v  =  kQ  z  cos  0.  A  change  in  u  will  effect  v  and  vice 
versa . 

From  this  point  on  it  becomes  clear  that  the  2-D 
case  is  not  simple  and  straightforward.  Equations  (A-32) 
and  (A-33)  are  not  independent  of  each  other.  Not  only 
will  there  be  a  variance  of  and  z^  but  also  a  covari¬ 
ance  . 

The  error  equation  for  the  one-dimensional  case 
described  a  two-dimensional  discriminator  curve  (i.e., 
for  an  error  z  a  A9  was  obtained) .  The  error  equations  in 
the  2-D  case  describe  a  three-dimensional  discriminator 
curve.  For  each  value  of  v  in  equation  (A-32) ,  there  is 
a  2-D  discriminator  curve  which  yields  a  Au  for  every 
Likewise,  in  equation  (A-33) ,  for  each  value  of  u  there  is 
a  2-D  curve  that  yields  a  Av  for  every  z^. 

The  variables  u  and  v  represent  the  actual  target 

A  A. 

position,  whereas  u  and  v  represent  an  estimate  of  the 
target  position.  For  computer  simulation  purposes,  the 
values  of  u  and  v  are  known  and  will  therefore  determine 
what  curve  in  (A-32)  and  (A-33)  to  use.  The  actual  error, 
however,  will  be  derived  from  equations  (A-30)  and  (A-31) 
and  the  resulting  and  z 2  are  applied  to  (A-32)  and 
(A-33)  to  obtain  Au  and  Av.  The  Au  and  Av  obtained  will 
theoretically  drive  the  boresite  to  the  target.  Because 
of  WGN,  the  boresite  will  not  point  directly  at  the  target. 
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Instead  it  will  be  within  some  variance  (or  covariance) 
of  the  actual  position.  A  new  error  will  develop  for  each 
observation  which  results  in  steering  the  boresite  closer 
to  the  target. 

The  input  to  the  azimuth  and  elevation  servo  must 
be  in  the  form  of  A 4>  and  A0  .  To  derive  these  values  note 
that 


and 


Au  =  k  x  sin  6  cos  d> 
o  T 


Av  =  kQ  z  cos  0 


(A-34 ) 


IA-35) 


If  we  let  k  x  =  kz-  1  we  have  two  equations  and  two 
o  o 

unknowns.  Au  and  Av  are  known  quantities  obtained  from 
the  discriminator  curves.  The  values  for  A0  and  A<J>  are 

A0  =  cor>"^  ( Av)  (A-36) 

Av  =  co,:  "^Au/sinfcos  (Av)))  (A-37) 

Since  0  and  4>  are  the  difference  between  boresite  and 
target,  0  =  A0  and  <p  -  A<£. 

With  this  information  a  tracking  loop  can  be 
implemented.  The  tracking  error  in  azimuth  and  elevation 
can  be  computed.  Recall  that  both  error  equations  contain 
azimuth  and  elevation  information.  It  should  also  be  noted 
that  the  results  for  A0  and  A(j>  are  dependent.  Besides  a 
separate  error  for  azimuth  and  elevation,  a  combined  error 
needs  to  be  developed. 
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This  appears  to  be  intuitively  consistent  with 
what  one  would  expect  from  a  two-dimensional  Amplitude- 


Comparison  monopulse.  In  Figure  A-4  the  2-D  beam  is  pic¬ 
tured  with  a  target  (x)  shown  off  boresite.  The  antenna 
functions  would  weight  such  a  target  return  in  both  azimuth 
and  elevation. 

Boresite 

EL 


Fig.  A-4.  A  Two-Dimensional  Beam  Showing  a 
Target  (x)  Off  Boresite 

In  order  to  move  the  boresite  quickly  and  most  efficiently, 
a  combined  azimuth-elevation  error  would  be  necessary. 

Steps  to  implement  a  two-dimensional  tracking  loop 
were  not  taken.  In  order  to  find  the  threshold  of  such  a 
system  and  quantify  it,  a  2-D  Z2,  and  CR  bound  is  necessary. 
This  is  the  obstacle  that  prevented  the  2-D  case  from  being 
pursued  further. 


The  Two-Dimensional  Lower  Bounds 

The  Cramer-Rao  (CR)  bound  is  defined  as  [Ref  6] 

Var(u)  >  [ - j - ~ - 1  (A-38 

E(iJ#  p  <r,  r  /A,u) }  . 

for  the  one-dimensional  case.  It  can  be  shown  that  for 
the  2-D  case  [Ref  6] 


Cov(u,v)  t* 


a2£n 

3u3v 


(A-39 


p(r1r2r3/A,u,v) } 


The  Ziv-Zakai  (ZZ)  bound  is  not  so  easily  trans¬ 
formed  into  two  dimensions.  The  one-dimensional  ZZ  bound 
is  stated  here  as  [Ref  3] 

e2(0)  max  sin2(0)P  (-0,0)  (A-40 

0<6<6  e 
- P 


where  0  is  defined  by 
P 


sin  0  = 

P 


sin  9  ,  +  sin  0 

max 


To  develop  a  two-dimensional  bound,  the  same  steps 
as  outlined  by  Ziv  and  Zakai  must  be  taken  [Ref  12] . 

The  first  obstacle  to  overcome  in  developing  a 
2-D  bound  is  finding  the  probability  of  error  P(e).  The 
P(e)  as  used  in  the  ZZ  bound  is  the  error  probability  of 
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the  best  procedure  for  deciding  whether  a  target  is  at  -6 
or  6  when  it  is  known  that  it  is  at  ore  of  these  positions 
and  each  has  equal  possibilities  [Ref  8] .  If  we  want  to 
expand  the  P(e)  to  two  dimensions,  four  regions  must  be 
considered. 

Instead  of  just  the  case  where  a  target  is  at  u^ 
or  Ujr  there  are  now  four  cases  shown  in  Figure  A-5. 
Instead  of  a  binary  detection  problem  we  now  have  an  M-ary 
detection  p^.'b _ a  where  M  =  4  [Ref  6j  . 


Fig.  A-5.  Decision  Regions  for  a 
Two-Dimensional  Case 


The  four  hypotheses  are 
H^s  r^  =  A/E  +  ni  (A-41) 

r2  =  A/I  ?2(U1V1)  +  n2 
r,  =  A/E  g_(u,v..)  +  n,. 


(A-4  2 ) 


H2:  rl  =  Ay/E  gl  ^U2V1^ 

r2  =  A/E  g2 (u2v1) 

r3  =  A/1  g3(u2Vl) 

H3:  r1  =  A/E  g1 (u2v2) 

r2  =  A/E  g2 (u2v2) 

r3  =  A/E  g3(u2v2) 

H4:  rl  =  Ay/E  gl(ulv2) 

r2  =  A/E  g2  (u-^v^ 

r3  =  A/E  g3(u1v2) 

From  detection  theory, 
costs  the  decision  region  R^  is 
for  which  [Ref  6:Ch.  5] 


+  n, 


+  n. 


+  n. 


+  n. 


(A-4 3) 


+  n. 


+  n. 


+  n. 


(A-44) 


+  n. 


+  n. 


with  probability  of  error 
given  as  the  values  of  r . 


P{H_.}p(r/A..)  >  P{Hk}p(r/Hk)  for  j#  k  (A-45) 

Equal  probability  of  hypotheses  is  assumed  so  equation 
(A-4  5)  becomes 

ptr/IK)  >  p(r/Hk)  (A-46 ) 

where  p(r/II.)  -  p  (r.  r_r_/H  , )  etc. 

3  1  2  3  j 
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For  decision  region  we  have 

p(r/H1)  >  p(r/H2) 
p(r/H1)  >  p(r/H3) 
p(r/H1)  >  p(r/H4) 

Similarly,  for  regions  R2,  R^ , 

R2:  p(r/H2)  >  p(r/H1) 
p(r/H2)  >  p(r/H3) 
p(r/H2)  >  p(r/H4) 

r3:  p(r/H3)  >  ptr/H^ 
p(r/H3)  >  p(r/H2) 
p(r/H3)  >  p (r/H4 ) 

r4:  p(r/H4)  >  p(r/H1) 
p(r/H4)  >  p(r/H2) 
p(r/H4)  >  p(R/H3) 


(A-4  7a) 
(A-47b) 
(A-4  7c ) 

and  R4 

(A-48a) 

(A-48b) 

(A-48c) 

(A-49a) 
(A-4 Ob) 
(A-49c) 

(A-50a) 

(A-50b) 

(A-50c) 


The  next  step  is  to  find  the  likelihood  ratios  for 
each  decision  region.  For  notational  purposes,  let 


P  (r/H  . ) 

Vr)  =  P(r/H.) 


S  p (r/H  . ) 

£jk(r)  =  plr/H^T 

Using  this  notation,  for  the  four  regions  we  have 

Rx:  £2<r)  <  0  R2:  £2(r)  >  0 

l}(r)  <  0  l23(r)  >  0 

i4iri  <  0  l24(r)  >  0 

R3:  £3 (r)  >  0  R4:  &4(r)  >  0 

£32(r)  >  0  &4 2 )  >  0 

H34(r)  >  0  J43(r)  >  0 

As  in  Chapter  II  of  this  thesis,  symmetry  will  be 
assumed  in  the  error.  The  probability  of  error  will  there¬ 
fore  be  the  same  in  each  decision  region.  The  total  P(e) 
will  be  the  sum  of  the  P(e) '  s  in  each  region.  This  result 
is  shown  as 

P  (e )  =  p{H1}P£(R1)  +  P{H2}P£(R2) 

+  P{H3)Pe(R3)  +  P{H4}P£(R4)  (A-51) 

Since  P{H..}  =  1/4  and  P£ (R^)  =  P£(R2)  =  pe(R3)  =  1e(R4) 
then 

P (e)  =  P  (R-i  )  (A-52) 

e  1 

The  statistics  for  the  likelihood  ratios  for  region  1  will 
be  determined  and  used  to  describe  the  P(e). 

-4 

A  ?• 
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By  following  the  same 


steps  outlined  in  Chapter  II 


we  get 


l  r 

:  t 
!  '« 
t  ; 


where 


and 


Also 


E{A2(r)/Hx}  =  ~  Eo(1-P2) 


(A-53) 


EQ  =  [g.12(u1v1)  +  912(U2V1)  +  g2  (ulVl) 

+  922(u2v1)  +  ^^iV  +  932(u2V1)3 
p2  =  ^-[gi(u1v1)g1(u2vi)  +  g2(uivi)92(u2vi> 

z  o 

+  g3(uivi)g3(u2vi)] 


EU2(r)/H2)  =  “  ^  e0(1“P2) 


(A-54) 


Var  (£  _  (r )  /H,  )  =  VarU2<r)/H  )  =  Nq/2  Eo(1~P2)  (A"55) 


a/e 


E{A3(r)/H1}  =  —  e0(1“p3) 


a/e 


E{ £^(r) /H3) }=  -  —  e0<1_P3^ 


(A-56) 

(A-57) 


Var  (£.,  (r)  /H,)  =  Var  U3  (r)  /H3>  =  NQ/2  EJI-P3)  (A-58) 


and 


a/e 


E{i4(r)/H1}  =  —  Eq(1-P4) 


a/e 


E{i4  (r)/H4)  - - 2~  e?o(1"P4) 


(A-59) 

(A-60) 
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A>  <<.  <wn>I 


Jbro 


VarU4(r)/H1)  =  Var  U4  (r)  /H4 )  =  N0/2Eq(1-p4)  (A-61) 


The  P  (e  )  is 


p(ji2>o/h2)p{h2}  +  pu3>o/h3)p{h3) 


+  P  (2,4>0/H4 )  P  {H4  ) 


(A-6  2) 


Since  P{H2>  =  P{H3>  =  P{H4>=  1/4#  equation  (A-62)  becomes 


P(e)  =  1/4P(£2>0/H2)  +  l/4P(i3>0/H3) 


+  1/4P(  >0/H4) 


Each  term  in  equation  (A-63)  can  be  expressed  as  a 
Q- function . 

EE  , 

P(e)  =  1/4{(1-Q[A(^j2(1-p2))s]) 

EE  , 

+  (l-Q[A(2^(l-p3))^]) 


(A-63) 


EE  , 

+  ( 1-Q  (A  ( 2jj~(!~  P4 )  )  )  > 

o 

Since  symmetry  is  assumed,  p2  =  P3  =  P4  =  P 


EE  , 

P  (  e)  =  3/4  {(l-Q[A(^(l-p))*]  )  } 

o 


(A-64 ) 


(A-6  5) 


Now  that  a  probability  of  error  has  been  deter¬ 
mined  the  next  step  is  to  derive  the  ZZ  bound.  Ziv  and 
Zakai  used  the  Tchebycheff  inequality  to  transform  their 
probability  of  error  into  a  form  including  a  mean-square 
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error.  Since  this  is  a  two-dimensional  case  a  two- 
dimensional  inequality  of  the  Tchebycheff  type  is  needed. 
Reference  [7]  shows  several  probability  inequalities  of 
the  Tchebycheff  type.  A  two-dimensional  one  called  the 
Berge  inequality  is  stated  as 


P  (max{ 


}  >  X) 


1  +y/  l-p2 


(A-70) 


where  and  x2  are  random  variables 

ui  =  E^} 

°i2  = 

°12  =  eUx^MXj-Uj)} 

p  =  — “ 

°1°2 

and  X  >  0 

As  can  be  seen  from  the  definition  of  p,  the  lower 
bound  on  the  covariance  could  possibly  be  obtained  from 
this  inequality;  however,  attempts  to  do  so  were  unsuccess¬ 
ful.  It  is  recommended  that  the  Berge  inequality  be 
investigated  further  as  to  its  possible  use  in  determining 
the  ZZ  bound  in  two  dimensions.  A  search  for  other  two- 
dimensional  inequalities  may  prove  worthwhile  also. 

Once  a  suitable  inequality  is  found  the  2-D  ZZ  bound 
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Appendix  C 

Plots  of  Moving  Target  MSE  Compared  to  the 
Cramer-Rao  and  Ziv-Zakai  Bounds 
from  35  dB  to  -10  dB 
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>  k  «  W*^>*i*  vW4» '  -?v*tt  fcW<-«i  t»  _■  r  \,  »--  4—.  ‘*:*Mi±*4ui&TtSaBiUt+M* '"ff< 


PROGRAM  MAIN< INPUT  *  OUTPUT ,L I  ST ,PL0T , T APE  5  =  1 \PUT , T APE6  =  L I  ST , 

1  TAPE9  =  PLCd 

C  *  THIS  PROGRAM  SIMULATES  A  STATIONARY  AND  MOV If  G  T  ARGET  « 

C  *  BEING  TRACKED  BY  AN  AMPLITUDE-COMPARISON  MONCPULSE  3  AD Ar , 

C  *  CALCULATES  THE  MEAN-SQUAPt  ERROR*  THE  CRAMER-PAC  BOUND* 

C  •  THE  ZIV-ZAKAI  BCUNO,  ANO  PLOTS  THE  RESULTS. 

COMMON/GRAFS/TITLE( 31 ,Y1T(21,Y2T(21 »Y3T(21,Y4T(?1 
COMMON /IN OU /KIN  *KOU  T  *  KP 

DIMENSION  TIME<52>  *  TPK(52I ,TR (52) *RMSEAC15*50 ) *R MS (15,50 ) 
DIMENSION  M100  1»SNRVC151 
REAL  NPWR*M*E2 
INTEGER  PULSE 
DOUBLE  PRECISION  DSEEO 
C  '*  INPUT  ALL  PARAMETERS 

DATA  KIN,KCUT*KP/5»6,9/ 

DATA  DEPRAD,PI*THEr AH/57.3,3.14159,0 ./ 

DATA  VT,R0/F09.N,69 OC ./ 

DATA  SNRV/35. ,3 0 . , 2 5 . ,20 . ,1 5. * 1 0 . ,7 . ,5 . , 3 . , 1 . , -2 . , -5 . ,-7 . , - 1 0 . , 

1  -2  0./ 

DATA  YIT/lCHMEAf.  SQUAR,10HE  ERROR  / 

DATA  Y3T/1GHCRAMER-PAG,10H  HOUND  / 

DATA  Y4T/1CHZIV-ZAKAI  ,10HBOUND  / 

REWIND  6 
REWIND  9 

CALL  PLOT S(0o*0.»9) 

5  CONTINUE 

WRITE*,-  MOVING  OR  STATIONARY  TARGETd  OR  21- 
READ*, IMS 

IF(IMS.LE.C)  GO  TC  510 

WRITE*,-  C-R  OR  Z-Z  BOUND  (1  0-  21- 

READ*,  IB 

WRITE*,-  ENTER  MSE(C-R!  ANO  Z-Z  BOUND  SCALES" 

READ*,TM,  TZ 

WRITE*,"  ENTER  DO  LOOP  -ANGE" 

READ*, II, 12 
IFdB.EQ.21  GO  TO  20 
Y2T (1 i=Y3T(l! 

Y2T (21=Y3T(21 
GO  TO  25 

20  Y2T(ll=Y4T(li 
Y2T (2  1  =  Y4T(Z1 
25  CONTINUE 

DSEEO=21047C0C0C.D0 
WRI TE(S,52! 

C  •  TOTAL  OBSERVATION  INTERVAL 
T=10Q.E-3 
TX=T/5.e-3 

C  *  SERVO  PARAMETERS 
SIG=600 . /TX 
0M=1483./TX 
NR  = 100 
T  RK  (  5 1 1 =0  . 
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TRK  <52)=TM/A. 

TR(  51  )  =  C. 

IF(IB.EQ.l)  TR(52)  =  TMM. 

I F( I B«  CQ*  2 )  TP<52)=TZ/4. 

T  IM  E( 5 1 ) =0  . 

TIME<52)=50 .+TX/6. 

DO  500  IP=I1»I2 
SNR  =  SNRV<  IP  ) 

SNR0  =  10.**<SNR/10.) 

SNRR  =  SNR0 
ROLD=R0 
SPUR=l.E-3 
DO  200  J= 1 » 1 5 
1=0 

IFCIMS.EQ.2)  GO  TO  30 

BRS ITE=1 • 5 

TRGT=0. 

GO  TO  35 
30  BRS I T  E=0 • 

TRG  T=l. 

35  CONTINUE 

THETA  =  TRGT-RFSITE 

C  *  CALL  GGNML  AND  SET  UP  THE  OBSERVATIONS  R1  AND  R2 
CALL  GGNML(DSEED,NR »P ) 

DO  100  PULSE=1 t 50 
PUL$EC=FL9AT(PULS£) 

RNE  U=SQRT(RI'  +  +  2  +  (VT*PlJL3EC*T  )*  *  2  ) 

IF( IHS.EQ.l )  SNRR=< (ROLD/RNEU)* *4 >*SNRC 

NPUR=SPWR/S\RR 

S=S  QRT <  NPUr  ) 

1=1  *1 
N1=R<I)*S 
1=1+1 
N2=R< I)*S 

GD=SIN($IN(THETAHn 
GS=COS(SIN<THETAH  )) 

THETA  =  THETA/OEPc  AD 
Rt  =  COS(SINC  THETA) )+Nl 
R2=  SI N( SI N< THETA)  >+\!2 
C  *  THE  ESTIMATOR  AND  THE  01 SCRIM  I  NAT  OR 
E=R  1+GD-R2+GS 
DTHCTA  =  -E*DEPc,AO 
IF< DTHETA.GT.1.5)  0THETA=1.5 
IF( OTHETA.LT,-l*5)  DTHETA=-1.5 
C  *  THE  SERVO 

SERVO=(l.-EXP(-SIG*  DM <SlG/CM)*SIN(OM*T)*COS COM*T  ))  ) 
OTHETP=DTHETA*SERVC 
C  *  MOVE  THE  ANTENNA 

BRS  ITE  =  0RSI  TE  +  DTr)ETP 
C  *  THE  EQUATION  TO  MOVE  THE  TARGET 
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I  Ft  IMS.EQ.l  )  TRGT  =  A  TANC  ( V T* PULSEC * T  )  /R 0  )  *  OE  PR  AD 
C  *  ADJUST  THETA  TO  REFLECT  NEW  POSITION  OF  TARGET  WITH  RESPECT 
C  *  TO  THE  BORESITE. 

THE  TA=TRGT-HR  SITE 
C  *  CALCULATE  MSE 
TRK  ER  =  -THE' A 
RMS  CA<J»PULSE>=TRKE°**2 
IFt THETA. GT. 1.5)  THtTA  =  1.5 
IFt THETA.LT. -1.51  THETA=-1.5 
THE TA=THETA/ DEPP  AD 
I  Ft  I  B.EQ*  2 )  GO  ’■C  A  0 
SM=1./SNRR/CCSCTHET A)**2 
GO  TO  45 
40  CONTINUE 

C  *  CALCULATE  CR  AND  ZZ  BCUNOS. 

UX=SINCVHETA> 

R0=C03<2.*UX> 

AX=SQRTt(  l.-RO/NPWR) 

YX=AX 

CALL  MONOR<YX*PX> 

PEX=1.-PX 
SM=PEX*UX**2 
45  RMS t Jf PULSE)=SM 

TIMECPULSE)-TX*PUL3EC 

theta=theta*oepfao 

SNR =1 0 « *ALCG10 ( SNPR  I  _ 

WRI TE(o»50)  BRSITE*TRGT *  PUL  SEC ♦TRKER»RMSEA(J ♦PULSE ) 

100  CONTINUE 
.  200  CONTINUE 

WRI TE<6«53> 

DO  400  K=l,50 
SUM=Q. 

SUM  1  =  0  « 

00  300  L= 1 » 15 
SUM  =SUM*R  MS  EA  (L  tK) 

SUM  1  =  SUM1 ♦RMS (L  *K  > 

300  CONTINUE 

TRK  CK)  =  SUM/15. 

TRC  K>  =  SUM1/15, 

WRI  TE  1 6  *  5  4)  TP.K  CK  >  *  SNR  *RC 
400  CONTINUE 

WRITE**"  ENTER  T I T L E ( 3 A 1 0 )• *  SNR  V 1 1 P ) 

READ(5»55>  TITLE 
C  *  PLOT  THE  RESULTS 

CALL  ORAWCCt  TIME»TRK»TR»Y3*Y4,2  *0*50* 1  ) 

500  CONTINUE 
GO  to  5 
510  CONTINUE 

WRITE**"  *  *  *  ROUTE  'PLOT*  TO  ON-LINE  PLOTTER" 

WRITE**"  *  *  *  ROUTE  •LIST*  TO  PRINTER" 

CALL  PLOT  E{  NDUM  ) 


on  o  n  non  on  on 


REWIND  6 
REWIND  9 

50  f  ORMATdOX,  M3X,E11 .5)  > 

52  F0RMAT(15X."BRS!TE"  » 9  X  »  " T  RG  T  "  •  9  X  ♦  "P  UL SEC"  »9  X  »  "  T  RKE  R "  ♦  ?  X  t  "  MSE  T«  > 

53  FQRMAT(16X*"MSEE"»1  L  X  »  "SNR"  »  12  X  ♦  HR'>  "  ) 

54  FORMAT(10X.3C3XtEU.5>> 

55  FORN*AT(3A10  ) 

56  F  OR  MAT ( 2 A1C  ) 

STOP 

END 

SUBROUTINE  DR  AWCC  C  X  » Y 1  »  Y2  ,  Y  3  »  Y4  ♦  NG»  A  PHS  t  NSC  ALE  ,  NP  :  I  NT  S  ,  JO  I  ) 

C OM  MON/GR  AF  S/ T I TLE C  3 ) »  Y1 TI T L E  <2 ) *  Y2  T I T LE  C 2 ) »  Y  3T I TLE ( 2 ) *  Y4 T I  TEE  <  2  ) 
COM  MON/ 1 N CU /KI N »KGU T ♦KPUNCH 

DIMENSION  XCl)*Yl<l>»Y2<l)*Y3<n»Y4(l),BIG<4>fSMALL<4>fXTITlEC4) 

- CHECK  CALL  PARAMETERS  FOR  ERRORS 

I  BO  OBQO=0 

IFCN GRAPHS «LE*4)G0  TO  7 
URITE(KOUTtlO) 

I  BO  OBOO  =  1 

7  IF( NSCALE.LE.NGF APHSJGO  TO  9 
URI TE(KOUT» 11  ) 

IBOCBOO=l 

9  XFCJOZN  .GE.  -1  .ANO.  JOIN  . LE »  1>GC  TO  14 
URI TE (KOUT » 13  ) 
lBOOBOO=l 

- ERROR  MESSAGES  FOP  CALL  PARAMETERS 

10  F0RMAT(1H-»1CX»*????7  TOO  MANY  ORDINATES  F CR  D^AUCC  CCCCCO 

11  FORMAT(lH-tlOX»*??7??  *  NSCALE*  LARGER  THAN  NGRAPHS  CCCCC*) 

13  FORMAT(1H-,10X. *?????  ’JOIN*  IS  OUTSIOE  ALLOWABLE  RANGE  CCCC*) 


14  CONTINUE 

- AUTOMATIC  SPECIFICATIONS 

0R0=4.  $  AB SC  =6 •  $  SIZE  =  .25  $  LETTERS  =  30 

NX=9$XTITLE (1 )  =  6H  PULSES 
NY1 =NY2=NY3=NY4-2C 


30  I  F(  I BOOBO  0  .NE.  1)GC  TO  33 
URI TE (KOUT »  27 ) 

27  FORMAT(lH-»5Xt*ERROR  IN  CALL  TO  DRAUCC  -  -  NO  PLOTS*//) 
RETURN 

- CENTER  THE  Y-AXES 

33  BOT  TOM=  <1 l.-ORD-S IZ  E  > /2 • 

C  -  —  -FIX  ORIGIN  FOR  Yl 

GRAPHS  =  NGRAPHS-1 
'  CALL  PLOT (GRAPHS»B0TT0M*-3) 

C 
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on  ">  no  oooo 


i  Ft  nscale  .eq.  c  >go  TC  34 

c  - CALCULATE  SCALE  FACTORS  FOR  EACH  ARRAY 

CALL  SCALE<X»ABSCt\PGINTS»l ) 

IF< NSCALE  • EE •  1 >G3  TO  80 
CALL  SCALE( Y1 ♦ORO,NPOINTS,l) 

IFt NGRAPHS. FQ.l >  GO  TO  34 
CALL  SCALE<  Y2  fOPDf  iPOINTStl  > 

IFt NGRAPHS. EQ. 2 >  GO  TO  34 

31  CALL  SCALEt Y3 ,0 PD, N PO I N  TS »1 ) 

IF( NGRAPHS. EQ. 3 >  GO  TC  34 

32  CALL  SCALEC Y4 ,GRD, N PCI  NTS , 1 > 

GO  TO  34 

- THIS  SECTION  ACCOMPLISHES  IOENTICAL  SCALING 

-—-SEARCH  Y1C  )  ANO  Y2<  )  FOR  MAX  AND  MIN  VALUES 
80  BIG (1 )=SMALLtl)=Yl C 1) 

BIG  t2)=SMALLt 2>=Y2t  1) 

DO  85  1=2, N POINTS 

IF(YKI)  •  GT  •  BIGtl)  )  BIG(1>=Y1(I> 

IF(Y2CI)  .GT.  BIG<2)  )  BIGC2)=Y2(I) 

IFCY1CI)  .LT.  SMALLC1)  )  SMALLC  1  )  =Y1  C  I  ) 

IF(Y2(I)  .LT.  SMALL (21  )  SMALL (2 ) =Y2 C I ) 

85  CONTINUE 

X  MA  X=  AM AX  1 1  BIG(1 ! • HI G ( 2 )  ) 

XMI  N=AMIN1(SMALL< 1) »SRALL(2)  ) 

I Ft NSCALE  .LT.  3>G0  TC  100 

-  —  -SEARCH  Y3<  )  FOR  ITS  MAX  AND  MIN  VALUES 
BIGt3)=SMALLt3)=Y3tl) 

DO  90  1=2  »h POINTS 

I  Ft  Y3  C I  I  .GT.  HI G ( 3  )  )  BIG<3)  =  Y3(I) 

IF(Y3CI)  .LT.  SMALL  C  3  >  )  SM ALL< 3 ) =Y3 C I ) 

90  CONTINUE 

IFCXMAX  .LT.  HI G(3 )  )XMAX=BIG(3> 

I  Ft  XM  IN  .GT.  SM ALL  t  3  )  )  XMI N  =  SMALL 1 3 ) 

I  Ft  NSCALE  .LT.  4)G0  TO  100 

- SEARCH  Y4t  )  FOP  ITS  MAX  AND  MIN  VALUES 

BlGt4)=SMALLt4)=Y4tl) 

DO  95  1  =  2  tNPOINTS 

I  Ft  Y4  1 1 )  .G^.  BI Gt4  )  )  BIG<4)=Y4(I> 

I  Ft  Y4  1 1 >  .LT.  SMALL  ( 4  )  )  SMALL  1 4 )  =  Y4  1 1 > 

95  CONTINUE 

I  Ft  XM AX  .LT.  BIGt4>  )XMAX  =  BIGt4) 

IFtXMAX  .LT.  BI  G 1 4  )  )XMAX=BIGt4) 

I  Ft  XM  IN  .GT.  SMALL  t  4  )  )  XM I N  =  SMA LL 1 4 ) 

-—-FIND  THE  ADJUSTED  MIN  AND  INCREMENT  PER  INCH 
100  BIGtl )=XMIN 
.  BIGt2)=XMAX 
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non  on  n  on  non 


CALL  3CALE(DIG»CPD.2  »1  > 

XMIN  =  BIG!3> 

DEL  T A=BIG!4  ) 

-—-PUT  THE  SCALING  INFORMATION  INTO  THE  ARRAYS 
GO  TO  (34»llb*ll0»105)tNSCALE 

105  Y4!  NPCINT$*1)=XMIN 
Y4t NPOINTS*?>=OELTA 
110  Y3( NPOINTS+l >=XMIN 
Y3C  NPOINTS+2)=.nELT  A 
115  Y2!  NPOINTS+1 > =XMIN 
Y2!NPOINTS*2)= DELTA 
Y 1 !  NPOINTS+1  >=X.MIN 
Yl! NPOI NTS+2 ) =DELT  A 

- IS  ANY  INDEPENDENT  SCALING  SEEDED? 

IFCNSCALE  .EQ.  2  .AND.  NGRAPHS  -GT.  2>G0  '0  31 
I F( NSC ALE  .EQ.  3  -AND.  NGRAPHS  .EQ.  4 ) GO  10  32 

34  A=0RD-.l 

-—-DRAM  AM  AXIS  AND  A  SYMBOL  FOR  THAT  AXIS 

CALL  AXIS<0.tO.,flTITLE*NYl  .  ORO  »SC  .  *  Y  1  (  NPCI  f.TS  +  1  )  »  Y  1  (  NPOI  NTS +  2  )  ) 

IF< JOIN.EQ.O )GC  TO  35 

CALL  SYMBOL!-. 4, A,. I5tl»90. .-1) 

35  IF! NGRAPHS. EQ.l >  GO  TO  51 

CALL  AXIS(-1.  ,0 .»Y2 TITLE tNY2*CRD. 90. tY2(NP0INTS  +  l )»Y2( NPOI NT $  +  2) ) 

IF!  JOIN.EC.OGO  TO  36 

CALL  SYMBCL!-1.4,A,  .15.2«SC  .  ,-1  > 

36  IF!  NGRAPHS'„C0.2>  GO  TO  51 

CALL  AXIS  <-2.  fO  ..Y3TITLC»NY3»CRD»90.  ,  Y3  !  NPOI  NT  S* 1 )  *  Y3  !  NP  01  N"  S*  2  )  ) 

IF! JOIN.EQ.O )GC  TC  37 

CALL  SYMBOL  !-2.4,  A*  .15,5,90  .«-l  ) 

37  IF! NGRAPHS. EQ. 3 )  GO  TO  51 

CALL  AXIS!-3.«C.fY4TITLEtNY4 » CP  0*90. *  Y 4  !  NPC I  NT S  «-l >  *  Y4 ! NPOI N ' S* 2 ) ) 

IF! JOIN.EQ.C ) GC  TC  51 

CALL  SYMBOL!-3.**A*.15*4t?G.,-l> 

51  CALL  AXISIO .  *0. * XT  I TLE » -NX » A BSC , 0 . » X ! NPOI NT$* 1 ) *  X ! NPO I N T S+2 ) ) 

- OECIOE  HOW  MANY  POINTS  SHOULD  BE  MARKED  UITH  A  SYMBOL  . 

IF! JOIN  .Nf.  01  GO  TO  72 

J=Q  . 

GO  TO  53 

72  IF! JOIN  .EQ.  1 )  GO  To  68 
J=-l 

WRI  TE  !  KOU  T  »  71  ) 

71  FORMAT! 1 H- » 5X »* EVER Y  POINT  IS  MAPKED  WITH  A  SYMBOL*) 

GO  TO  53 
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no  o  n  o  no 


69  NUM=ABSC*.25  ♦  1. 

IFINUM  .LT.  3>NUM=3 
J=NPOINTS/NUM 
IFtJ  .LT.  1  >0  =  1 

69  JJ=IAB3(J) 

WRI  TE <  KOUT  t 70 )  JJ 

70  F0RMAT<1H-»5X**EVERY* *I4t*  TH  POINT  IS  HARKED  WITH  A  SYMBOL* ) 
- PLOT  X  VS  Y3 

53  CALL  LINE(X»Yl*t'.POI,\TS»l*d*l  ) 

IFCNGRAPHS  .EQ.  1>G0  TO  60 
- PLOT  X  VS  Y2 

56  CALL  LINE(X.Y2*  *.P0  INTStl«J«2) 

I F( NGR  AP  HS  .EQ.  2  >  G  0  TO  60 
- PLOT  x  VS  Y 3 

57  CALL  LINEIX»Y3»  f.POI  NTS»1»J»5) 

I FC  NGRAPHS  .EQ.  3>GC  TG  60 
- PLOT  X  VS  Y A 

58  CALL  LINECX *Y4tNP0INTS»l »J.4> 

- SEE  IF  A  TITLE  IS  DESIRED 

60  IFC LETTERS. EQ.O )GO  TO  61 

C  CENTER  TITLE  OVER  GRAPH 

XX=0. 

CALL  SYMB0L<XX«0RD*.425»SIZE»TITLE»G..LETTERS> 

C  -  —  -ORAW  BOX  AROUND  PLOT 

RLM=-<GRAPHS*1.  ) 

RRM  =.ABSC+i  . 

RRM=-1. 

RTM=ORD+l  . 

CALL  PLOT  (RLM»RBM»3  ) 

CALL  PLOT  (R  R  M*RBM  *2  ) 

CALL  PL0T(RRM,RTM,2  ) 

CALL  PLOT (RLM»RTM»2  ) 

CALL  PLOT  (SLMtRBM*2  ) 

C  - POSITION  PEN  FOR  THE  NEXT  CALL  DRAWCC 

61  CALL  PLOT (AHSC+5. »-BOTTOM»-3  ) 

RETURN 

END 
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